clear; clc; close all;

%% 参数设置
a = 0.1;
b = 0.1;
k = 1.85;
c_values = [0, 0.2];  % 分别绘制两种情况

% 定义 q0 的变化范围
q0_min = -20;
q0_max = 20;
num_q0 = 500;             % 定义 num_q0
q0_list = linspace(q0_min, q0_max, num_q0);
x0_fixed = 0.1;           % 固定初始 x

% 迭代参数
transient = 1000;         % 舍去暂态步数
N_plot = 100;             % 用于绘制分岔图的步数
N_LE = 1000;              % 用于计算 Lyapunov 指数的有效迭代步数

%% 对不同的 c 分别计算并绘图
for ci = 1:length(c_values)
    c = c_values(ci);
    
    % 初始化保存数据的变量
    bifurcation = cell(1, num_q0);
    LE1_vals = zeros(1, num_q0);
    LE2_vals = zeros(1, num_q0);
    
    % 对每个不同的 q0 初始条件进行迭代
    for idx = 1:num_q0
        q0 = q0_list(idx);
        % 初始化状态
        x = x0_fixed;
        q = q0;
        % 先迭代 transient 步，去除暂态
        for i = 1:transient
            [x, q] = mptm(x, q, a, b, c, k);
        end
        % 记录后续 N_plot 步的 x 值作为分岔图数据
        x_temp = zeros(1, N_plot);
        for i = 1:N_plot
            [x, q] = mptm(x, q, a, b, c, k);
            x_temp(i) = x;
        end
        bifurcation{idx} = x_temp;
        
        % 计算 Lyapunov 指数
        [LE1, LE2] = LEs(@mptm, x0_fixed, q0, a, b, c, k, N_LE, transient);
        LE1_vals(idx) = LE1;
        LE2_vals(idx) = LE2;
    end
    
    %% 绘制分岔图和 Lyapunov 指数
    figure;
    subplot(2,1,1);
    hold on;
    % 绘制每个 q0 对应的 x 值（分岔图）
    for idx = 1:num_q0
        q_val = q0_list(idx);
        x_data = bifurcation{idx};
        plot(q_val*ones(size(x_data)), x_data, '.k', 'MarkerSize', 1);
    end
    xlabel('q0');
    ylabel('x');
    title(sprintf('Bifurcation Diagram (c = %g)', c));
    grid on;
    
    subplot(2,1,2);
    plot(q0_list, LE1_vals, 'b', 'LineWidth', 1.5);
    hold on;
    plot(q0_list, LE2_vals, 'r', 'LineWidth', 1.5);
    xlabel('q0');
    ylabel('Lyapunov Exponents');
    legend('LE1','LE2','Location','Best');
    title(sprintf('Lyapunov Exponents (c = %g)', c));
    grid on;
end
